Last updated: 2016-08-30
Code version: fe93f23076b12cabdcee6348f3a79dca62c1de88
We have collected dendritic cells (DCs) from two populations: the first are individuals that have latent TB infections (putatively resistant), and the second are individuals that have recovered from an active TB infection (putatively sensitive). We performed RNA-seq on MTB-infected DCs and mock-treated controls. Because of our small sample size, we are not performing an eQTL analysis. Instead, we are performing a differential expression (DE) analysis to find gene expression differences between susceptible and resistant individuals.
To test for DE, we used the following linear model (implemented with limma+vooom):
\[ Y\ \sim\beta_0+X_{treat}\beta_{treat}+X_{status}\beta_{status}+X_{treat,status}\beta_{treat,status}+I+\epsilon \]
where \(Y\) is the expression level of a gene, \(\beta_{treat}\) is the fixed effect of treatment with MTB, \(\beta_{status}\) is the fixed effect of susceptibility status, \(\beta_{treat,status}\) is the fixed effect of interaction between treatment and susceptibility status, and \(I\) is the random effect of individual (implemented via duplicateCorrelation).
The assumption is that genes which are differentially expressed in our in vitro system will have regulatory variants that affect TB susceptibility. In other words, we expect that SNPs nearby the differentially expressed genes will be enriched for low p-values obtained from GWAS of TB susceptibility. We would like to test this by combining the DE results with published GWAS results. Specifically we use the following framework.
If there is a decrease in the GWAS p-values for increasing DE effect sizes, this suggests that these genes contain nearby genetic variants that affects TB susceptibility.
In this analysis, we use the p-values from a GWAS of TB susceptibility in Ghana published in Thye et al., 2010. The initial results look promising, i.e. there is a decrease in GWAS p-values for increasing DE effect sizes. Ideally we would also test for an enrichment from additional TB GWAS.
library("limma")
library("data.table")
library("biomaRt")
library("SNPlocs.Hsapiens.dbSNP144.GRCh38")
library("dplyr")
library("GenomicRanges")
library("ggplot2")
library("cowplot")
Input the results of the model fit by limma.
fit <- readRDS("../data/results-limma-fit.rds")
Input summary statistics from Thye et al., 2010.
gwas_thye_ghana <- fread("../data/OUT_PLINK_Ghana.txt", data.table = FALSE,
verbose = FALSE)
##
Read 4.1% of 2436894 rows
Read 39.4% of 2436894 rows
Read 50.5% of 2436894 rows
Read 87.0% of 2436894 rows
Read 2436894 rows and 6 (of 6) columns from 0.079 GB file in 00:00:07
The file contains the chromosome, SNP rsID, odds ratio (OR), and p-value (PVAL). Because the file does not contain any allele frequency data, it should not be possible to identify study participants (Craig et al., 2011).
colnames(gwas_thye_ghana)
## [1] "CHR" "SNP" "TEST" "NMISS" "OR" "PVAL"
Assign a GWAS summary statistic to each gene
11336 genes were tested for differential expression.
gene_names <- rownames(fit$coefficients)
head(gene_names)
## [1] "ENSG00000279457" "ENSG00000188976" "ENSG00000187961" "ENSG00000187583"
## [5] "ENSG00000187642" "ENSG00000188290"
Obtain the transcription start site (TSS) for each gene.
# Ensembl 83, Dec 2015, grch38.p5, hg38
# ensembl <- useMart(host = "dec2015.archive.ensembl.org",
# biomart = "ENSEMBL_MART_ENSEMBL",
# dataset = "hsapiens_gene_ensembl")
tss_all_fname <- "../data/tss-all.rds"
if (file.exists(tss_all_fname)) {
tss_all <- readRDS(tss_all_fname)
} else {
tss_all <- getBM(attributes = c("ensembl_gene_id", "chromosome_name",
"transcription_start_site", "strand"),
filters = "ensembl_gene_id",
values = gene_names,
mart = ensembl)
saveRDS(tss_all, file = tss_all_fname)
}
head(tss_all)
## ensembl_gene_id chromosome_name transcription_start_site strand
## 1 ENSG00000000419 20 50958550 -1
## 2 ENSG00000000419 20 50958555 -1
## 3 ENSG00000000419 20 50945861 -1
## 4 ENSG00000000419 20 50958521 -1
## 5 ENSG00000000419 20 50958532 -1
## 6 ENSG00000000457 1 169893952 -1
This returns the TSS for each transcript of gene. To reduce this to one number, take the most upstream TSS (i.e. the most 5’ for genes on positive strand and the most 3’ for genes on the negative strand).
tss <- tss_all %>%
group_by(ensembl_gene_id) %>%
summarize(chr = chromosome_name[1],
strand = strand[1],
tss = if (strand == 1) min(transcription_start_site) else max(transcription_start_site))
head(tss)
## Source: local data frame [6 x 4]
##
## ensembl_gene_id chr strand tss
## (chr) (chr) (int) (int)
## 1 ENSG00000000419 20 -1 50958555
## 2 ENSG00000000457 1 -1 169894267
## 3 ENSG00000000460 1 1 169662007
## 4 ENSG00000000938 1 -1 27635277
## 5 ENSG00000000971 1 1 196651878
## 6 ENSG00000001036 6 -1 143511690
The window to search for SNPs for a gene will be +/- 50 kb from the TSS.
window <- 50000
tss$start <- tss$tss - window
tss$end <- tss$tss + window
tss$strand <- ifelse(tss$strand == 1, "+", "-")
tss_gr <- makeGRangesFromDataFrame(tss, keep.extra.columns = TRUE)
seqlevels(tss_gr) <- paste0("ch", seqlevels(tss_gr))
tss_gr
## GRanges object with 11336 ranges and 2 metadata columns:
## seqnames ranges strand | ensembl_gene_id
## <Rle> <IRanges> <Rle> | <character>
## [1] ch20 [ 50908555, 51008555] - | ENSG00000000419
## [2] ch1 [169844267, 169944267] - | ENSG00000000457
## [3] ch1 [169612007, 169712007] + | ENSG00000000460
## [4] ch1 [ 27585277, 27685277] - | ENSG00000000938
## [5] ch1 [196601878, 196701878] + | ENSG00000000971
## ... ... ... ... ... ...
## [11332] ch16 [ 29765952, 29865952] + | ENSG00000280789
## [11333] ch7 [ 30500761, 30600761] - | ENSG00000281039
## [11334] ch10 [ 13560047, 13660047] + | ENSG00000282246
## [11335] ch1 [111453760, 111553760] - | ENSG00000282608
## [11336] ch7 [150350702, 150450702] + | ENSG00000283013
## tss
## <integer>
## [1] 50958555
## [2] 169894267
## [3] 169662007
## [4] 27635277
## [5] 196651878
## ... ...
## [11332] 29815952
## [11333] 30550761
## [11334] 13610047
## [11335] 111503760
## [11336] 150400702
## -------
## seqinfo: 25 sequences from an unspecified genome; no seqlengths
Convert rsID to genomic coordinates.
snp_coords_fname <- "../data/snp-coords.rds"
if (file.exists(snp_coords_fname)) {
snp_coords <- readRDS(snp_coords_fname)
} else {
snp_coords <- snpsById(SNPlocs.Hsapiens.dbSNP144.GRCh38, gwas_thye_ghana$SNP,
ifnotfound = "drop")
saveRDS(snp_coords, file = snp_coords_fname)
}
stopifnot(mcols(snp_coords)$RefSNP_id %in% gwas_thye_ghana$SNP)
Overlap the SNPs with the genes.
overlaps <- findOverlaps(snp_coords, tss_gr, ignore.strand = TRUE)
How many SNPs were found for each gene?
snps_per_gene <- countSubjectHits(overlaps)
stopifnot(length(snps_per_gene) == length(gene_names))
summary(snps_per_gene)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 24.00 52.00 55.77 81.00 268.00
sum(snps_per_gene > 0)
## [1] 10268
How many genes were found for each SNP?
genes_per_snp <- countQueryHits(overlaps)
stopifnot(length(genes_per_snp) == length(snp_coords))
table(genes_per_snp)
## genes_per_snp
## 0 1 2 3 4 5 6 7 8
## 1976866 315960 91371 25511 7465 2589 1005 523 270
## 9 10 11 12 13
## 162 56 12 11 4
Convert overlap results to use original SNP and gene names.
results <- data.frame(as.matrix(overlaps))
colnames(results) <- c("rsID", "gene")
results$rsID <- mcols(snp_coords)$RefSNP_id[results$rsID]
results$gene <- mcols(tss_gr)$ensembl_gene_id[results$gene]
head(results)
## rsID gene
## 1 rs6603811 ENSG00000215790
## 2 rs6603811 ENSG00000008130
## 3 rs7531583 ENSG00000215790
## 4 rs7531583 ENSG00000008130
## 5 rs12044597 ENSG00000215790
## 6 rs12044597 ENSG00000008130
Add GWAS p-values.
rownames(gwas_thye_ghana) <- gwas_thye_ghana$SNP
results$gwas_p <- gwas_thye_ghana[results$rsID, "PVAL"]
stopifnot(!is.na(results$gwas_p))
For each gene, assign the minimum p-value of all its nearby SNPs.
results <- results %>%
group_by(gene) %>%
summarize(gwas_p = min(gwas_p),
n_snps = n())
As expected, there is a negative correlation between the minimum GWAS p-value and the number of SNPs. However, there is no technical reason that genes assigned many SNPs should be more likely to be differentially expressed in our in vitro infection of dendritic cells with MTB.
cor(results$n_snps, results$gwas_p)
## [1] -0.4070354
plot(results$n_snps, results$gwas_p, xlab = "Number of SNPs nearby gene",
ylab = "Minimum GWAS p-value",
main = "Relationship between number of tested SNPs near gene\nand the minimum GWAS p-value of these SNPs")
Add coefficients from differential expression analysis.
limma_coef <- fit$t
results <- merge(results, limma_coef, by.x = "gene", by.y = "row.names")
Convert the coefficients to their absolute value.
for (test in colnames(limma_coef)) {
results[, test] <- abs(results[, test])
}
Create categorical variables to identify genes with absolute effect sizes greater or less than 1.
for (test in colnames(limma_coef)) {
new_col <- paste0("de_", test)
results[, new_col] <- ifelse(results[, test] > 1, "|logFC| > 1",
"|logFC| <= 1")
results[, new_col] <- factor(results[, new_col],
levels = c("|logFC| <= 1", "|logFC| > 1"))
}
results %>% select(starts_with("de_")) %>% summary
## de_diff_before de_diff_after de_treat_resist
## |logFC| <= 1:5556 |logFC| <= 1:6777 |logFC| <= 1: 865
## |logFC| > 1 :4712 |logFC| > 1 :3491 |logFC| > 1 :9403
## de_treat_suscept de_diff_treat
## |logFC| <= 1:1514 |logFC| <= 1:6276
## |logFC| > 1 :8754 |logFC| > 1 :3992
As a sanity check to confirm that our strategy of choosing the minimum SNP per gene is not biased, check that the mean expression level (i.e. the intercept term from the model) is not associated with lower p-values.
mean_expression_level <- fit$Amean[results$gene]
stopifnot(results$gene == names(mean_expression_level))
results$mean_expression_level <- mean_expression_level
test_intercept <- lm(gwas_p ~ mean_expression_level, data = results)
summary(test_intercept)
##
## Call:
## lm(formula = gwas_p ~ mean_expression_level, data = results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.10550 -0.06898 -0.04377 0.01520 0.90465
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.071624 0.002940 24.360 < 2e-16 ***
## mean_expression_level 0.002808 0.000582 4.824 1.43e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1267 on 10266 degrees of freedom
## Multiple R-squared: 0.002262, Adjusted R-squared: 0.002165
## F-statistic: 23.27 on 1 and 10266 DF, p-value: 1.426e-06
This is confirmed. In fact, reassuringly it is the exact opposite pattern: genes with higher mean expression (and thus with more statistical power to call differential expression) are associated with higher GWAS p-values.
Is there an enrichment of GWAS p-values for SNPs nearby genes that are differentially expressed following treatment with MTB in resistant individuals?
dist_treat_resist <- ggplot(results, aes(x = treat_resist, y = gwas_p)) +
geom_point(alpha = 0.5) +
geom_vline(xintercept = 1, color = "red") +
labs(x = "|logFC| treat_resist", y = "GWAS p-value",
title = "GWAS p-value vs. DE effect size")
wilcox_treat_resist <- wilcox.test(gwas_p ~ de_treat_resist, data = results)
enrich_treat_resist <- ggplot(results, aes(x = de_treat_resist, y = gwas_p)) +
geom_boxplot() +
labs(title = sprintf("Wilcox test p-value: %.2e", wilcox_treat_resist$p.value),
x = "DE effect size", y = "GWAS p-value")
plot_grid(dist_treat_resist, enrich_treat_resist, nrow = 1)
Is there an enrichment of GWAS p-values for SNPs nearby genes that are differentially expressed following treatment with MTB in susceptible individuals?
dist_treat_suscept <- dist_treat_resist %+% aes(x = treat_suscept) +
labs(x = "|logFC| treat_suscept")
wilcox_treat_suscept <- wilcox.test(gwas_p ~ de_treat_suscept, data = results)
enrich_treat_suscept <- enrich_treat_resist %+% aes(x = de_treat_suscept) +
labs(title = sprintf("Wilcox test p-value: %.2e", wilcox_treat_suscept$p.value),
x = "DE effect size")
plot_grid(dist_treat_suscept, enrich_treat_suscept, nrow = 1)
Is there an enrichment of GWAS p-values for SNPs nearby genes that are differentially expressed between individuals that are resistant or susceptible to TB in the noninfected state?
dist_diff_before <- dist_treat_resist %+% aes(x = diff_before) +
labs(x = "|logFC| diff_before")
wilcox_diff_before <- wilcox.test(gwas_p ~ de_diff_before, data = results)
enrich_diff_before <- enrich_treat_resist %+% aes(x = de_diff_before) +
labs(title = sprintf("Wilcox test p-value: %.2e", wilcox_diff_before$p.value),
x = "DE effect size")
plot_grid(dist_diff_before, enrich_diff_before, nrow = 1)
Is there an enrichment of GWAS p-values for SNPs nearby genes that are differentially expressed between individuals that are resistant or susceptible to TB in the infected state?
dist_diff_after <- dist_treat_resist %+% aes(x = diff_after) +
labs(x = "|logFC| diff_after")
wilcox_diff_after <- wilcox.test(gwas_p ~ de_diff_after, data = results)
enrich_diff_after <- enrich_treat_resist %+% aes(x = de_diff_after) +
labs(title = sprintf("Wilcox test p-value: %.2e", wilcox_diff_after$p.value),
x = "DE effect size")
plot_grid(dist_diff_after, enrich_diff_after, nrow = 1)
Is there an enrichment of GWAS p-values for SNPs nearby genes in which the differential expression following treatment with MTB is different between individuals that are resistant or susceptible to TB?
dist_diff_treat <- dist_treat_resist %+% aes(x = diff_treat) +
labs(x = "|logFC| diff_treat")
wilcox_diff_treat <- wilcox.test(gwas_p ~ de_diff_treat, data = results)
enrich_diff_treat <- enrich_treat_resist %+% aes(x = de_diff_treat) +
labs(title = sprintf("Wilcox test p-value: %.2e", wilcox_diff_treat$p.value),
x = "DE effect size")
plot_grid(dist_diff_treat, enrich_diff_treat, nrow = 1)
Is there an increase in GWAS SNPs nearby genes that are differentially expressed following treatment with MTB?
dist_snps_treat_resist <- ggplot(results, aes(x = treat_resist, y = n_snps)) +
geom_point(alpha = 0.5) +
geom_vline(xintercept = 1, color = "red") +
labs(x = "|logFC| treat_resist", y = "Number of SNPs nearby gene",
title = "Number of SNPs vs. DE effect size")
wilcox_snps_treat_resist <- wilcox.test(n_snps ~ de_treat_resist, data = results)
enrich_treat_resist <- ggplot(results, aes(x = de_treat_resist, y = n_snps)) +
geom_boxplot() +
labs(title = sprintf("Wilcox test p-value: %.2e", wilcox_snps_treat_resist$p.value),
x = "DE effect size", y = "Number of SNPs nearby gene")
plot_grid(dist_snps_treat_resist, enrich_treat_resist, nrow = 1)
Is there an increase in GWAS SNPs nearby genes that are differentially expressed between individuals that are resistant or susceptible to TB in the noninfected state?
dist_snps_diff_before <- dist_snps_treat_resist %+% aes(x = diff_before) +
labs(x = "|logFC| diff_before")
wilcox_snps_diff_before <- wilcox.test(n_snps ~ de_diff_before, data = results)
enrich_diff_before <- enrich_treat_resist %+% aes(x = de_diff_before) +
labs(title = sprintf("Wilcox test p-value: %.2e", wilcox_snps_diff_before$p.value),
x = "DE effect size")
plot_grid(dist_snps_diff_before, enrich_diff_before, nrow = 1)
Is there an increase in GWAS SNPs nearby genes in which the differential expression following treatment with MTB is different between individuals that are resistant or susceptible to TB?
dist_snps_diff_treat <- dist_snps_treat_resist %+% aes(x = diff_treat) +
labs(x = "|logFC| diff_treat")
wilcox_snps_diff_treat <- wilcox.test(n_snps ~ de_diff_treat, data = results)
enrich_diff_treat <- enrich_treat_resist %+% aes(x = de_diff_treat) +
labs(title = sprintf("Wilcox test p-value: %.2e", wilcox_snps_diff_treat$p.value),
x = "DE effect size")
plot_grid(dist_snps_diff_treat, enrich_diff_treat, nrow = 1)
I obtained significant results when splitting the data in two: genes with effect size greater than and below 1. If I perform a linear regression, only the effect of treatment is statistically significant. It is easier to see why this is when viewing the log-log transformed plot, which is always the plot in the right panel below.
Effect of treatment in resistant individuals.
summary(lm(gwas_p ~ treat_resist, data = results))
##
## Call:
## lm(formula = gwas_p ~ treat_resist, data = results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.08856 -0.06998 -0.04371 0.01534 0.90219
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0886181 0.0019924 44.478 < 2e-16 ***
## treat_resist -0.0005235 0.0001953 -2.681 0.00735 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1268 on 10266 degrees of freedom
## Multiple R-squared: 0.0006996, Adjusted R-squared: 0.0006023
## F-statistic: 7.187 on 1 and 10266 DF, p-value: 0.007355
plot_grid(dist_treat_resist + geom_smooth(method = "lm"),
dist_treat_resist + geom_smooth(method = "lm") + scale_x_log10() + scale_y_log10())
Effect of treatment in susceptible individuals.
summary(lm(gwas_p ~ treat_suscept, data = results))
##
## Call:
## lm(formula = gwas_p ~ treat_suscept, data = results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.08729 -0.06998 -0.04393 0.01532 0.90194
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0873487 0.0020099 43.458 <2e-16 ***
## treat_suscept -0.0006633 0.0003613 -1.836 0.0664 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1268 on 10266 degrees of freedom
## Multiple R-squared: 0.0003282, Adjusted R-squared: 0.0002308
## F-statistic: 3.37 on 1 and 10266 DF, p-value: 0.06642
plot_grid(dist_treat_suscept + geom_smooth(method = "lm"),
dist_treat_suscept + geom_smooth(method = "lm") + scale_x_log10() + scale_y_log10())
Effect of susceptibility status in noninfected state.
summary(lm(gwas_p ~ diff_before, data = results))
##
## Call:
## lm(formula = gwas_p ~ diff_before, data = results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.08446 -0.07003 -0.04407 0.01531 0.90224
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.447e-02 2.109e-03 40.057 <2e-16 ***
## diff_before -9.263e-06 1.606e-03 -0.006 0.995
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1269 on 10266 degrees of freedom
## Multiple R-squared: 3.24e-09, Adjusted R-squared: -9.741e-05
## F-statistic: 3.326e-05 on 1 and 10266 DF, p-value: 0.9954
plot_grid(dist_diff_before + geom_smooth(method = "lm"),
dist_diff_before + geom_smooth(method = "lm") + scale_x_log10() + scale_y_log10())
Effect of susceptibility status in infected state.
summary(lm(gwas_p ~ diff_after, data = results))
##
## Call:
## lm(formula = gwas_p ~ diff_after, data = results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.08877 -0.06961 -0.04366 0.01468 0.90691
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.088904 0.002038 43.616 < 2e-16 ***
## diff_after -0.005261 0.001906 -2.761 0.00578 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1268 on 10266 degrees of freedom
## Multiple R-squared: 0.0007419, Adjusted R-squared: 0.0006446
## F-statistic: 7.622 on 1 and 10266 DF, p-value: 0.005776
plot_grid(dist_diff_after + geom_smooth(method = "lm"),
dist_diff_after + geom_smooth(method = "lm") + scale_x_log10() + scale_y_log10())
Effect of interaction term.
summary(lm(gwas_p ~ diff_treat, data = results))
##
## Call:
## lm(formula = gwas_p ~ diff_treat, data = results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.08549 -0.07000 -0.04405 0.01556 0.90299
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.085625 0.002133 40.149 <2e-16 ***
## diff_treat -0.001269 0.001883 -0.674 0.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1269 on 10266 degrees of freedom
## Multiple R-squared: 4.423e-05, Adjusted R-squared: -5.317e-05
## F-statistic: 0.4541 on 1 and 10266 DF, p-value: 0.5004
plot_grid(dist_diff_treat + geom_smooth(method = "lm"),
dist_diff_treat + geom_smooth(method = "lm") + scale_x_log10() + scale_y_log10())
sessionInfo()
## R version 3.2.1 Patched (2015-07-12 r68650)
## Platform: x86_64-unknown-linux-gnu (64-bit)
## Running under: Scientific Linux release 6.7 (Carbon)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] cowplot_0.6.2
## [2] ggplot2_2.1.0
## [3] dplyr_0.4.3
## [4] SNPlocs.Hsapiens.dbSNP144.GRCh38_0.99.11
## [5] BSgenome_1.36.3
## [6] rtracklayer_1.30.4
## [7] Biostrings_2.38.4
## [8] XVector_0.10.0
## [9] GenomicRanges_1.22.4
## [10] GenomeInfoDb_1.6.3
## [11] IRanges_2.4.8
## [12] S4Vectors_0.8.11
## [13] BiocGenerics_0.16.1
## [14] biomaRt_2.26.1
## [15] data.table_1.9.6
## [16] limma_3.26.9
## [17] knitr_1.13.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.5 plyr_1.8.3
## [3] formatR_1.4 futile.logger_1.4.1
## [5] bitops_1.0-6 futile.options_1.0.0
## [7] tools_3.2.1 zlibbioc_1.14.0
## [9] digest_0.6.9 gtable_0.2.0
## [11] RSQLite_1.0.0 evaluate_0.9
## [13] DBI_0.4-1 yaml_2.1.13
## [15] stringr_1.0.0 grid_3.2.1
## [17] Biobase_2.28.0 R6_2.1.2
## [19] AnnotationDbi_1.32.3 XML_3.98-1.4
## [21] BiocParallel_1.2.22 rmarkdown_1.0
## [23] lambda.r_1.1.7 magrittr_1.5
## [25] scales_0.4.0 Rsamtools_1.22.0
## [27] htmltools_0.3.5 GenomicAlignments_1.6.3
## [29] assertthat_0.1 SummarizedExperiment_1.0.2
## [31] colorspace_1.2-6 labeling_0.3
## [33] stringi_1.1.1 lazyeval_0.1.10
## [35] munsell_0.4.3 RCurl_1.95-4.8
## [37] chron_2.3-47